#include <iostream>
#include <fstream>
#include <cmath>
#include <cstdio>
#include <sstream>

#include "shekel_function.h"
#include "griewank_function.h"
#include "rosenbrock_function.h"
#include "zakharov_function.h"

#include "Time.h"
#include "HybridMonteCarlo.h"
#include "ZFunction.h"
#include "wrapperFunctions.h"


#define PI   3.1415926535897932385

using namespace std;

void run_HMC(int argc, char *argv[ ], Function *func, string& directory)
{
  char dummy2[80];
  char dummy3[80];
  char dummy4[80];
  int seed;
  double TempStart, TempStep, TempEnd;
  double SuccessRatio;
  int    StepNumber; // Nmd
  int	   NSample;	  // Nmc number of samples for each run
  double StepLength; // = 0.01; //1e-2;
  int NEachTempRun; // number of runs for each temperature
  int N; // number of degrees of freedom
  int nMin;     // every nMin trajectory is to be stored in the container
  //double *mass;




  /////////////////////////////////////////////////////
  // checks to see if 11 arguments have been specified
  // if not then reads them from a file
  /////////////////////////////////////////////////////

  if (argc == 15)
  {
    // required content of input arguments
    //
    // argv[1]     seed
    // argv[2]     temperature start
    // argv[3]     temperature step
    // argv[4]     temperature end
    // argv[5]     metropolis success ratio
    // argv[6]     N_md    (StepNumber)
    // argv[7]     N_mc    (NSample)
    // argv[8]    start step length (StepLength)
    // argv[9]    N_runs  (NEachTempRun)
    // argv[10]    the number of dimensions

    seed = atoi(argv[1]);
    TempStart = atof(argv[2]);
    TempStep = atof(argv[3]);
    TempEnd = atof(argv[4]);
    SuccessRatio = atof(argv[5]);
    StepNumber = atoi(argv[6]);
    NSample = atoi(argv[7]);	 
    StepLength = atof(argv[8]); 
    NEachTempRun = atoi(argv[9]);
    N = atoi(argv[10]);
    strcpy_s(dummy2, 80, argv[11]);
    strcpy_s(dummy3, 80, argv[12]);
    nMin = atoi(argv[13]);
    strcpy_s(dummy4, 80, argv[14]);
  }
  else if (argc == 1) // if no arguments given to executable do this
  {
    ifstream file("input/params_hmc.txt");

    if (! file)
    {
      cout << "Cannot open params.in\n";
      exit(1);
    }

    string dummy;
    string dummy1;
    string dummy2;

    file >> seed >> TempStart >> TempStep
      >> TempEnd >> SuccessRatio
      >> StepNumber >> NSample >> StepLength
      >> NEachTempRun >> N >> dummy >> dummy1 >> nMin >> dummy2;

    //newA(mass,N);

    /*for(int i = 1; i<=N; i++)
    {
      file >> mass[i];
    }*/

    file.close();
  }
  else
  {
    cout << "Incorrect input.\n";
    exit(1);
  }



  ///////////////////////////////////////
  // save params used for HMC experiment
  // used for debugging. And to keep a 
  // record of the job executed
  ///////////////////////////////////////

  ofstream fpParams("./output/HMC_params.used");
  fpParams << "Seed               " << seed << endl;
  fpParams << "Temperature start  " << TempStart << endl;
  fpParams << "Temperature step   " << TempStep << endl;
  fpParams << "Temperature end    " << TempEnd << endl;
  fpParams << "metropolis         " << SuccessRatio << endl;
  fpParams << "N_md               " << StepNumber << endl;
  fpParams << "N_mc               " << NSample << endl;
  fpParams << "StepLength         " << StepLength << endl;
  fpParams << "N_runs             " << NEachTempRun << endl;
  fpParams << "N                  " << N << endl;
  fpParams << "nMin               " << nMin << endl;

  fpParams.close();

  initiate_ran1(seed);    // initiates the random number


  Time timer; // to measure calculation times
  cout << "prepare timer " << timer;


  double *Input;                // -> initial solution
  double *solLocalMinimized;    // -> solution obtained after local minimization
  double *solCont;              // -> local solution container



  ////////////////////////////////////////////////////////////////
  // Open files HMC_Solution.out and HMC_initial_values.out,
  // and prepare them for output. HMC_Solution.out holds the
  // best function value after each run both without local
  // minimisation and with local minimisation. HMC_initial_values
  // holds the initial coordinates (debugging).
  ////////////////////////////////////////////////////////////////

  ofstream outEndHMC_Solution( string(directory + "HMC_Solution.out").c_str());

  outEndHMC_Solution << "Dimension: " << N << endl
    << "Smallest function value after each run\n"
	  << "(LM = Local Minimization)\n\n"
	  << "     without LM           with LM           "/*Best in Cont           */
    << "Function Evaluations          Distance"
    << "          Run Number\n";

  ofstream outIniVal("output/HMC_initial_values.out");

  Function *MyFunction;       // Creates an instance of the function
  MyFunction = func;          // Sets the function = to the function passed as an argument

  newA(Input, N);
  newA(solLocalMinimized,N);
  newA(solCont, N);

  
  
  //ofstream rejected("output/nRejected.reject");



  ////////////////////////////////////////////
  // creates a HMC object with the specified
  // arguments; get the time; and begin a run
  ////////////////////////////////////////////

  HybridMonteCarlo MySampler(N, StepNumber, StepLength, MyFunction, 1.0 - SuccessRatio);
  //HybridMonteCarlo MySampler(N, StepNumber, StepLength, MyFunction, 1.0 - SuccessRatio, mass);

  //deleteA(mass);

  double fret;              // used for local CGM minimisation
  int Nrun = 1;             // which run

  cout << "Time: " << timer; 

  // If the start-temperature specified is different from the finish-temperature
  // then do a HMC run for each of the temperature from TempStart toTempEnd
  // in steps of TempStep
  for (double MainTemp = TempStart; MainTemp <= TempEnd+1e-07; MainTemp+=TempStep)
  { 
    double gloBest = 1e16;
	  MySampler.SetTemp(MainTemp);        // Sets the temperature in the HMC object
    outEndHMC_Solution << "\nTemperature = " << MainTemp << endl << endl;

    // For each temperature you may do multiple runs if the
    // input parameter corresponding to NEachTempRun > 1
	  for(int i=1; i<=NEachTempRun; i++)  
	  {
      // Output to screen the temperature and the run number
      // for that temp
      cout << "T = " << MainTemp << "; i = " << i << endl;

      MyFunction->nCall = 0;
      // Populate the Input array which a random start position/configuration
		  MyFunction->createRandomSol(Input);

      // Print out the temperature, run number and initial
      // starting configuration (coordinates)
		  outIniVal << "\nTemperature = " << MainTemp << " : Run No. " << i << endl << endl;
      outIniVal << "Dimension | Starting values" << endl;
		  for (int k=1; k<=N; k++)
      {
        outIniVal << setw(5) << k << setw(18) << Input[k] << endl;
      }

      // reset tracker before starting a new HMC run
		  MySampler.ResetTracker(Input);             // calculates the function value of sol

		  // used to diagnose the program (only used in CreateSamplesDiagnosis)
		  /*char funcSample[80]; 
		  sprintf_s(funcSample, 80, string(directory + "HMC_diagnosis%d.out").c_str(), Nrun);
		  char coordSample[80]; 
		  sprintf_s(coordSample, 80, string(directory + "HMC_coord%d.coord").c_str(), Nrun);*/

      MySampler.CreateSamples(Input, NSample, nMin);
      //MySampler.CreateSamplesDiagnosis(Input, NSample, funcSample, coordSample); // do a run consisting of NSamples
                              // (i.e. trajectories) each of StepNumber 
                              // Molecular dynamics steps

      // print out best solution (before local minimisation)
		  /*char resultHMC[80];

		  sprintf_s(resultHMC, 40, string(directory + "result_HMC%d.txt").c_str(), Nrun);
      MyFunction->print_sol(Input,ChiSq,resultHMC,0,0);*/

		  //int lNAccept, lNTotal; // lNAccept=number of accepted MC moves. lNTotal=total number of MC moves
		  //MySampler.GetAcceptAndTotal(lNAccept, lNTotal);
		  //cout << "\nSave files: " << resultCTS << ", ";

      // print out best solution so far (without local minimisation)
		  double lBestFuncValue = MyFunction->BestSoFar(solLocalMinimized);
      double lBestInCont = 1e16;

		  outEndHMC_Solution << setw(14) << lBestFuncValue;

      wrapperCGM(solLocalMinimized, N, 1e-3, &lBestFuncValue, MyFunction, directory);


      // Used for MRSA
      char buf [4];

      if(nMin != 0)
      {
        vector<Sol>& cont = MyFunction->Result.getCoordFromCont();
        ofstream outCont( string(directory + "HMC_Container" + itoa(Nrun, buf, 10) + ".out").c_str() ); 
        
        outCont << "Nmc     Func Val Bef CGM         Func Val Aft CGM" << endl;

        for (unsigned int j = 0; j < cont.size(); j++)
        {
          outCont << MyFunction->Result.whichNmc[j] << setw(18) << cont[j].funcVal << setw(25);

          for (int i = 0; i<N; i++)
          {
            solCont[i+1] = cont[j].coord[i];
          }

          wrapperCGM(solCont,N,1e-3,&cont[j].funcVal,MyFunction,directory);  // calls local minimisation

          for (int m = 0; m<N; m++)
          {
            cont[j].coord[m] = solCont[m+1];
          }

          outCont << cont[j].funcVal << endl;

        }

        outCont.close();

        /*cout << "cont[0].coord[0] = " << cont[0].coord[0] << endl;
        cout << "cont[0].funcVal = " << cont[0].funcVal << endl;
        cout << "-----------------------------\n" << endl;*/

        MyFunction->Result.BestInCont(lBestInCont, solCont);

		    // print out best solution (after local minimisation)
		    //sprintf_s(resultHMC, 80, string(directory + "HMC_result_CGM%d.txt").c_str(), Nrun);
		    //MyFunction->print_sol(resultHMC);

        // after local minimisation lBestFuncValue has been changed
        // this sets lBestFuncValue equal to fret which is
        // then outputted
        fret = lBestFuncValue;
		    outEndHMC_Solution << setw(20) << fret << setw(20) << lBestInCont << setw(21) 
          << MyFunction->nCall << setw(26) << MyFunction->distanceFromGlobalMinimum(solCont) 
          << setw(16) << Nrun << endl;

        if(lBestInCont < fret)
        {
          if(gloBest > lBestInCont)
          {
            gloBest = lBestInCont;
          }
        }
        else if(lBestInCont > fret)
        {
          if(gloBest > fret)
          {
            gloBest = fret;
          }
        }
        else if(lBestInCont == fret)
        {
          if(gloBest > lBestInCont && gloBest > fret)
          {
            gloBest = fret;
          }
        }
      }
      else
      {
        fret = lBestFuncValue;
		    outEndHMC_Solution << setw(20) << fret << setw(21) 
          << MyFunction->nCall << setw(26) << MyFunction->distanceFromGlobalMinimum(solLocalMinimized) 
          << setw(16) << Nrun << endl;

        if(gloBest > fret)
        {
          gloBest = fret;
        } 
      }

		  // reset the step length (i.e. delta time used in the leap-frog algorithm)
		  // to the value originally specified as the 10 input parameter. This then 
      // implies that when executing multiple runs, each run start with the same 
      // value for the input parameter.
		  MySampler.NewStepLength(StepLength); 

      //rejected << Nrun << setw(13) << MySampler.nReject << endl;

		  Nrun++; // next run
    }

    outEndHMC_Solution << "\nGlobal best value over all runs:\t" << setw(9) << gloBest << endl << endl
      << "---------------------------------------------------------------" << endl;

  }

  outEndHMC_Solution.close();
  outIniVal.close();
  //rejected.close();


  // print final output to screen
  cout << endl << endl
       << "Output files have been written to " << directory << endl;
  cout << "The total time for the simulation was " << timer << endl << endl;
  cout << "---------------------------" << endl << endl;

}